Performance benchmarks:

As can be seen in the performance benchmarks, shifting from stdlib random to numpy.random.default_rng comes with a performance hit. I dug into this a bit more. The two main random calls that are made in these benchmark models are shuffle and choice. np.random.default_rng().shuffle() is about 20% faster than random.shuffle(). However, np.random.default_rng().choice is several orders of magnitude slower than random.choice(). In fact, for the specific case of choosing a single item randomly from a collection, it is advised to use collection[rng.integers(0, len(collection)] instead, but this is still twice as expensive as random.choice.

More broadly, numpy.random.default_rng() is designed for generating and operating on numpy arrays. So, calls to create a single random number (e.g., rng.random, rng.integers) are all much slower than their stdlib random alternatives. However, once the size becomes larger, they become increasingly faster than a list expression with stdlib random. So, shifting to rng comes with additional implementation considerations to work around this performance hit. See some of the updated examples for what that can look like.

Thoughts are welcome. In your opinion, is the improved quality of the random number generators and the fact that you have local state with numpy worth the performance hit?

Given how common it is to draw a single random number (i.e., random.random(), I was wondering about adding a shorthand for this and see how fast this was.

def test_rand(rng, initial_size=100):
    initial_values = rng.random(size=initial_size)

    while True:
        for entry in initial_values:
            yield entry
        initial_value = rng.random(size=initial_size)

def draw():
    return next(a)

a = test_rand(rng, 200)  

And then I timed these as shown below

%%timeit
[random.random() for _ in range(250)]
8.88 μs ± 72.4 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

%%timeit
[draw() for _ in range(250)]
20.7 μs ± 130 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

%%timeit
[next(a) for _ in range(250)]
16.4 μs ± 121 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

%%timeit
[rng.random() for _ in range(250)]
67.9 μs ± 273 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

So, a naive shift from random.random to rng.random, is almost an order of magnitude slower. We can effectively halve this by using the generator and advancing this generator directly via next, wrapping this next call up in another function adds some overhead, but it is still a lot faster than a naive replacement. So, I am considering adding a method to the model, model.rand, that performs this task. The other option would be to figure out how to subclass numpy.random.Generator and add this rand method to it, because then you can just do self.rng.rand as a much faster alternative for rng.random().

Performance benchmarks:

More broadly, numpy.random.default_rng() is designed for generating and operating on numpy arrays.

I don't know if it's a good idea, but: Can you pre-generate an array (like a 100 or a 1000), use those numbers, and everytime you used up the array generate a new one? Basically generate random numbers in batch?